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Abstract. We present a one-dimensional model of the formation and viscous evolution of protoplanetary disks. 
The formation of the early disk is modeled as the result of the gravitational collapse of an isothermal molecular 
cloud. The disk's viscous evolution is integrated according to two parameterizations of turbulence: The classical 
q representation and a parameterization, representative of non-linear turbulence driven by the keplerian shear. 
We apply the model to DM Tau and GM Aur, two classical T-Tauri stars with relatively well-characterized disks, 
retrieving the evolution of their surface density with time. We perform a systematic Monte-Carlo exploration of 
the parameter space (i.e. values of the a-/3 parameters, and of the temperature and rotation rate in the molecular 
cloud) to find the values that are compatible with the observed disk surface density distribution, star and disk 
mass, age and present accretion rate. We find that the observations for DM Tau require 0.001 < a < 0.1 or 
2 x 10~ 5 < 13 < 5 x 10~ 4 . For GM Aur, we find that the turbulent viscosity is such that 4 x 10" 4 < a < 0.01 
or 2 x 10 -6 < P < 8 x 10 -5 . These relatively large values show that an efficient turbulent diffusion mechanism is 
present at distances larger than ~ 10 AU. This is to be compared to studies of the variations of accretion rates 
of T-Tauri stars versus age that mostly probe the inner disks, but also yield values of a ~ 0.01. We show that 
the mechanism responsible for turbulent diffusion at large orbital distances most probably cannot be convection 
because of its suppression at low optical depths. 

Key words. Accretion, accretion disks, Solar system: formation, planetary systems, planetary systems: formation, 
planetary systems: protoplanetary disks 



1. Introduction 

The presence of circumstellar disks has been proposed long 
ago as a necessary preliminary of star formation, simply 
because the decrease of the moment of inertia of a collaps- 
ing molecular cloud core by a factor 10 16 prevents the di- 
rect formation of a central compact object (the star) with- 
out a significant loss of angular momentum. These disks 
are now routinely discovered, and both statistical infor- 
mation of disk properties obtained from the spectral en- 
ergy distributions and characterization of individual disks 
from direct measurements can be obtained (see Protostars 
& Planets IF for detailed reviews on the subject). 

The basic principle yielding the formation of a disk 
around a protostar also governs its evolution: material is 
allowed to be accreted from the disk onto the protostar 
only by decreasing its specific angular momentum. The 
conservation of mass and total angular momentum implies 
that some of the material in the disk must be transported 
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outward in order for accretion to be possible (e.g. Lynden- 
Bell and Pringle, 1974; Pringle, 1981; Cassen, 1994). 

One of the major puzzles in understanding circumstel- 
lar disks remains the mechanism by which angular mo- 
mentum is transported. It is easy to show that viscosity 
on a microscopic scale is unable to explain the dissipa- 
tion of disks on any reasonable timescale. Beyond that, 
the problem resides not in finding a mechanism of angular 
momentum transport but rather on deciding which one is 
appropriate and what are the magnitudes of the resulting 
angular momentum transport and associated heat dissipa- 
tion (see e.g. Cassen 1994; Lin and Papaloizou 1996; Stone 
et al. 2000). It is generally believed that disks of relatively 
small masses transport their angular momentum by tur- 
bulence and that their evolution can be calculated using 
an effective turbulent viscosity. The parameterization of 
this viscosity then allows the calculation of the evolution 
of circumstellar disks to be performed over the timescales 
of interest (Myrs). 

In this paper, we confront a 1-dimcnsional theoreti- 
cal model of the formation and evolution of circumstel- 
lar disks to direct observations of the disks around the 
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T-Tauri stars DM Tauri and GM Aurigae. These disks 
have the particularity that their outer surface densities 
(between <~ 50 and 800 AU) are constrained by millimeter- 
wavelengths observations (Guilloteau and Dutrey 1998; 
Dutrey et al. 1998), continuum IR emission from the dust 
(Kitamura et al. 2002), and observations of star light scat- 
tered by the disk (Schneider et al. 2003). We use two differ- 
ent parameterizations of the turbulent viscosity (named a 
and P, see §3 hereafter) and two models of molecular cloud 
collapse. 

The simplicity of our model allows for an extensive ex- 
ploration of the space of parameters. In total more than 
10,000 models were run. We can therefore constrain the 
magnitude of the turbulent viscosity necessary to produce 
disks of the observed age, mass, and surface density. The 
observational constraints however are subject to very large 
uncertainties and instead of one simple scenario of forma- 
tion and evolution for each disk we find several kind of so- 
lutions to the current state of DM Tau and GM Aur. Our 
results can be used to address the questions of the source 
of the turbulence, the initial conditions and the type of 
collapse of the molecular cloud. They are also helpful in 
providing a framework for planet formation, as all the im- 
portant quantities such as densities and temperatures can 
be calculated as a function of time and distance to the 
central star. 

The structure of this article is as follows: In Section 
2 we present the observational characteristics of the DM 
Tau and GM Aur systems, and in particular the param- 
eters that must be fitted by the models. The numerical 
model itself and the physics behind it (including the possi- 
ble sources of angular momentum transport) are presented 
in Section 3. Section 4 presents some model examples and 
the strategy employed to fit the observations. We present 
our selected sets of models in section 5 and discuss the 
consequences for the disks global evolution, the turbulent 
angular momentum transport and the initial type of col- 
lapse. Finally, section 6 summarizes the main results of 
the article. 

2. Observations of circumstellar disks: DM Tau 
and GM Aur 

DM Tau and GM Aur are both located in the Taurus- 
Auriga region, a young stellar association relatively close 
to Earth, <~ 140 pc (Kenyon et al. 1994), with a low abun- 
dance of massive stars and a relatively cold (~ 10 — 20 K) 
environment (van Dishoeck et al. 1993). These stars have 
been chosen for this study because their disks have been 
characterized at millimeter to ultraviolet wavelengths, and 
because the emission of CO and its isotomers has been 
detected by millimctric interferometers to large distances 
(500 to 800 AU) to the stars. For these disks, we thus 
have access to observations sensitive to the presence of 
dust (spectral energy distribution in the infrared and at 
millimetric wavelength, visible images of reflected stellar 
light, images in the millimetric continuum), to the pres- 
ence of hot inner gas accreted by the star (excess emission 



mainly in the U band) and to the presence of cold gaseous 
CO in the outer parts of the disks (interferometric 12 CO 
and 13 CO emission maps). 

2.1. DM Tau 

DM Tau is a low mass (0.50 M Q ) relatively old T-Tauri 
star (~ 5 x 10 6 yr). Its gas disk was discovered by 
Guilloteau and Dutrey (1994) with the IRAM 30-m tele- 
scope and independently by Handa et al. (1995). Further 
observations and a x 2 analysis (Guilloteau and Dutrey 
1998) allowed the determination of the disk inclination 
and showed a rotation rate consistent with a kcplerian 
profile. More recently, Dartois et al. (2003) combined ob- 
servations of CO lines of different C-isotopes to better in- 
fer the surface densities and also estimate the disk vertical 
temperature profile. Their study confirmed theoretical ex- 
pectations that the disk mid-plane at ~ 100 AU is colder 
(T - 15 K) than 2-3 scale heights higher (T - 30 K). The 
disk radius, 850 AU, and the slope of the density profile, 
p — 1.5 derived by Guilloteau and Dutrey differ from a 
previous study by Saito et al (1995) who derived 350 AU 
and p = 2.0 ± 0.3. The discrepancy can probably be at- 
tributed to the lower resolution (5" compared to 3") and 
sensitivity of the observations by Saito et al. 

The outer disk radius measured for DM Tau cither 
from the mm continuum emission (Kitamura et al. 2002) 
or from HST near IR coronographic images (Grady et al. 
2002) is assumed to be a lower limit because of the system- 
atically smaller lower sensitivities (in terms of equivalent 
surface densities) of these techniques compared to the CO 
observations. 

Recent observations of DM Tau in the near- to mid- 
infrared (2.9— 13.5 /im) and preliminary radiative transfer 
models by Bergin et al. (2004) seem to indicate that the 
inner region of the disk consists of three parts: (i) close 
to the star, and up to ~ 4AU, a flat, optically thin disk, 
cleared of most of its dust (and probably gas); (ii) At 
4AU, a "wall", heated by direct stellar irradiation to <~ 
120 K, and extending vertically to ±0.5 AU relative to the 
midplanc; (iii) A "standard" optically thick disk beyond 
that distance. 

2.2. GM Aurigae 

The evidence of the presence of a disk around GM Aur was 
first obtained from radiometric observations by Koerner 
et al. (1993). Higher resolution interferometric maps were 
later obtained by Dutrey et al. (1998), and Kitamura 
et al. (2002) using the same techniques as for DM Tau. 
Depending on the inclination of the disk, two solutions 
can be found for the star mass, i.e. M» = 0.5 or 0.9 M 
with most probable values close to 0.9. Its disk appears 
to be slightly smaller than that of DM Tau, and have 
an outer radius ~ 525 AU, as determined from 12 CO ob- 
servations. High-resolution coronographic HST images of 
scattered near-infrared light (Wood et al. 2002; Schneider 
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System: 




TlA/T Tan 
U 1V1 ± ail 


/~i A T Any. 


Distance: 




140 pc ± 10 % 


140 pc ± 10 % 


Age: 




1.5 - 7 Myr 


1 - 10 Myr 


Spectral Type: 




Ml 


K7 


Star Mass: 




0.4-0.6 M Q 


0.5-1.0 M Q 


Luminosity: 




0.25 Lq 


0.74 L 


Temperature: 




3720 K 


4060 K 


Star Radius: 




0.0065 AU 


0.0085 AU 


Constraints from CO emission lines: 


Rout 




850 ± 20 AU 


525 ± 20 AU 


^(lOOAU) 


> 


0.05 gcm~ 2 




S (flo„t) 


> 


1.6xl0~ 3 gcm~ 2 


3.0xl0~ 3 gcm~ 2 


Constraints from continuum dust 


emission: 


Rout 


> 


500 AU 


280 AU 


S(100Ai7) 


< 


4.3 gcm~ 2 


23.0 gem" 2 


£(100Ai7) 


> 


0.23 gem" 2 


1.4 gcm~ 2 


P 




0.48 - 1.03 


1.11 - 1.34 


Accretion rate: 






Log M 




-7.95 ±0.54 


-8.02 ±0.54 



[Meyr- 1 ] 



Table 1. Physical parameters for DM Tau and GM Aur: 
Stellar parameters are from Simon et al (2001), CO constraints 
from Dartois et al. (2003) and Dutrey et al (1998), dust emis- 
sion constraints from Kitamura et al. (2002) (with an addi- 
tional 1/3 to 3 uncertainty factor -see text) and accretion rates 
from Hartmann et al. (1998) and Gullbring et al. (1998). 



et al. 2003) indicate that the disk is outwardly flared and 
extends at least to 300 AU. 

As for DM Tau, the near- to mid-infrared observa- 
tions indicate the presence of an inner disk clearing within 
6.5 AU, with a 120 K "wall" at that distance and an opti- 
cally thick disk beyond (Bcrgin et al. 2004). The presence 
of such a gap (but a slightly smaller, i.e. ~ 4 AU one) was 
also an outcome of the modeling of a less accurate SED 
(Chiang and Goldrcich 1999, Wood et al. 2002), and at- 
tributed to the presence of a ~ 2 Mj planet orbiting at 2.5 
AU (Wood et al, 2002; Rice et al. 2003). This explanation 
remains mostly speculative at this point, however. 

2.3. Inferred constraints 

We detail hereafter the observational constraints that are 
directly relevant to the disk evolution calculations and 
that are summarized in Table 1. We choose to adopt 
relatively conservative constraints. (More restrictive con- 
straints may be obtained from SED fitting and direct anal- 
ysis of the emission maps which is beyond the scope of the 
present article.) 

Stellar masses: They are inferred from the keplerian ro- 
tation rates measured in 12 CO (Simon et al. 2001). 
The uncertainty on this parameter is essentially due 
to that of the distance given by Hipparcos to an accu- 



racy ~ 10%. The inferred stellar masses depend also 
on the the disk inclination. 

Ages: This is a crucial, but unfortunately very uncer- 
tain parameter entering the models. The uncertainty 
on the stellar masses propagates into that on the in- 
ferred age. Taking that into account, and using evolu- 
tion tracks from Baraffe et al. (1998, 2002), we derive 
for both stars a range of ages that covers ages pub- 
lished by different groups (e.g. Guilloteau & Dutrey 
1998; Hartmann et al. 1998). 

Accretion rates: As discussed by Hartmann et al. (1998), 
an accretion luminosity L acc can be derived from the 
measured excess emission observed in the U band. It 
is then turned into an accretion rate following the re- 
lation: 



where i?* and M* are the radius and mass of the star, 
respectively, and 7 is a factor that accounts for the 
distance from which the material free-falls onto the 
star. This parameter is expected to be <~ 0.8 due to 
the presence of a magnetic cavity inside 5 R+. The un- 
certainties on M stem both from those linked to the 
determination of the bolometric accretion luminosity, 
and those related to the uncertain 7. Hartmann et al. 
estimate the uncertainty on M to be within a factor 
3, with more optimistic error estimates leading to a 
factor ~ 2, given the likely presence of the magnetic 
cavity. 

Surface densities from measured CO emission maps: In 
the case of DM Tau, measurements of 12 CO, 13 CO 
and C 18 lines allow for a direct determination 
of the density profile of gaseous CO from orbital 
distances <~ 100 AU to the outer disk (Dartois et al. 
2003). Interestingly, the outer radius of the disk is 
different depending on the molecular species that 
is considered, which is interpreted as selective pho- 
todissociation in the outer disk. 12 CO being the most 
abundant isotomcr, it is also the most resistant to 
photodissociation and yields the largest radius (used 
for this work). Because CO furthermore tends to 
condense in these cold regions (e.g. Aikawa et al. 
1999), the gas mass can only be thought as a lower 
limit to the outer disk mass. A comparison of CO 
and dust emission in DM Tau indeed yields a disk 
mass that is smaller by a factor ~ 5 when inferred 
from the CO measurement than when inferred from 
the dust (Dartois et al. 2003). This appears to be 
consistent with CO condensation, but may also be 
partially explained by a larger continuum opacity. 
The same conclusions hold for GM Aur, but the set 
of observations is more limited, because it is based 
on 13 CO J = 1 -» (Dutrey et al. 1996) and 12 CO 
J = 2 -> 1 (Dutrey et al. 1998) emission. In Table 1, 
the minimum values of the surface density required to 
explain the CO emission are derived from Dartois et 
al. (2003) for DM Tau and Dutrey et al. (1998) for 
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GM Aur. Conservatively, we assumed that all CO was 
in gaseous form and with a solar abundance to derive 
these minimum densities. 
Surface densities from continuum emission maps: 

Because continuum millimetric and sub-millimetric 
emissions are optically thin, their measurements 
inform us on the surface density and global masses 
of the disks. However, for a given flux, the surface 
density (or equivalently disk mass) derived is inversely 
proportionnal to an uncertain opacity coefficient k v . 
Following Bcckwith et al. (1990), this opacity is often 
approximated by: 

where Ko ~ 0.1cm 2 g _1 and [3 are parameters that de- 
pend on the grain composition and size (e.g. Beckwith 
et al. 2000). ko is the opacity per total mass (gas+dust) 
and therefore depends on an assumed dust to gas ra- 
tio. The [3 parameter can be obtained from the slope of 
the spectral energy distribution in the mm region (e.g. 
Beckwith et al. 1990; Chiang & Goldreich et al. 1997). 
Uncertainties on k are often neglected in the calcula- 
tion of disk masses and surface densities although they 
are probably large (Beckwith et al. 1990; Pollack et al. 
1994; Hartmann et al. 1998; D'Alessio et al. 2001). 
In this work, we adopt the matches to the SEDs and 2 
mm emission maps obtained by Kitamura et al. (2002) 
to constrain the surface density at 100 AU. We choose 
not to use their constraints on the disk masses and 
radii and on the slope of the surface density profile be- 
cause the loss of sensitivity of the observations at large 
orbital distances and the limited spatial resolution im- 
ply that these are very model-dependant. Kitamura et 
al. do include a range of allowed (3 values to estimate 
the characteristics of the disks, but do not allow for 
variations in «o- In pratice, the values of [3 for DM Tau 
are smaller (by ~ 0.45) than those for GM Aur. Given 
the standard opacity law (eq. 1), this implies that the 
assumed 1.3 mm opacity is about twice larger for DM 
Tau than for GM Aur. On the other hand, models of 
grain growth indicate that low values of $ < 1 are 
consistent with the presence of large grains, and also 
generally lead to smaller opacities in the millimetric 
(Miyake & Nakagawa 1993; D'Alessio et al. 2001). This 
implies that the larger disk masses derived for GM Aur 
may be (at least partly) an artefact of the opacity law. 
In order to account for the uncertainty on kq we hence 
consider cases for which the 100 AU gas surface density 
is either multiplied or divided by a factor 3 compared 
to the values obtained by Kitamura et al. (2002). 

3. Model description 

Our aim in this paper is to obtain evolutionary models 
able to satisfy the observational constraints provided by 
Table 1 and perform a sensitivity analysis of the model pa- 
rameters for DM Tau and GM Aur. This requires solving 



the evolution of the gas disk for ^10 Myrs and search- 
ing not only for the "best" models but for the ensemble 
of parameters compatible with the observations and their 
uncertainties. With current (and near-future) computing 
facilities we must restrict ourselves to a relatively sim- 
ple system of ID radial equations in which all quantities 
have been vertically averaged and radiative transport is 
approximated. Consequences of this averaging appear to 
be relatively mild in terms of accuracy (Hure and Galliano, 
2001; see also Ruden and Lin, 1986). This is certainly not 
a concerning problem, given that we are seeking order- 
of-magnitudc constraints on the turbulent viscosity and 
other relevant parameters. 

We hereafter describe our system of equations in the 
framework of a simple, but relatively complete model in- 
cluding -hopefully- most relevant physical mechanisms. 
We also discuss its limitations in Section 3.5. 

3.1. Viscous evolution of the surface density 

The evolution of a disk follows conservation of angular 
momentum and mass. As mass is accreted into the inner 
region of the disk and finally onto the central star, part 
of the disk material must migrate outwards to conserve 
angular momentum. If the disk is axisymmetric and geo- 
metrically thin, the dynamics are only radially dependent 
and an effective turbulent viscosity, v, can be postulated 
as the angular momentum transport mechanism. In this 
case, the equation giving the evolution of the disk surface 
density, E, can be written as, 




where we have introduced S(r,t), a source term account- 
ing for the accretion from the molecular cloud core onto 
the circumstellar disk (see §3.3). 

The general solution of (2) is a disk that progressively 
accretes mass inwards while redistributing a part of the 
outermost material further away to conserve angular mo- 
mentum until all of the material is accreted and all of the 
angular momentum has been transported away (Lynden- 
Bcll and Pringlc, 1974). This process is modulated by the 
value of the (turbulent) viscosity v. 

3.2. Two parameterizations of turbulence: a and (3 
models 

A prescription for the turbulent viscosity was originally 
postulated by Shakura and Sunyaev (1973) and has been 
extensively used in models of turbulent disks known since 
then as alpha disk models. According to this prescription, 

v = ac s H. (3) 

The free parameter a controls the amount of turbulence 
in a turbulent medium where the scale height H, and the 
isothermal sound speed c s , are upper limits to the mixing 
length and turbulent velocity, respectively. Accretion rates 
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inferred in T Tauri stars are compatible with a w 10 -2 
(Hartmann et al, 1998). As a caveat, one should note that 
it is difficult to imagine any physical process that would 
yield a viscosity obeying strictly to this parameterization, 
with a being spatially uniform and temporally constant. 
In this sense, any value of a that can be retrieved has to 
be considered as an undefined average over both time and 
space of the adimensional quantity v/(c s H). 

A possible source of angular momentum transport that 
can be approximately parameterized in the a-framework is 
due to the so-called magnetorotational instability (Balbus 
and Hawlcy 1991, 1998, 2000). The instability arises from 
the fact that ionized elements in a magnetized environ- 
ment tend to conserve their velocity. In a keplerian mag- 
netized disk, an element being displaced inward will have 
a lower velocity than other elements at the same location. 
Therefore, it will also have a smaller angular momentum 
and it will tend to sink further in (sec Tcrqucm 2002 for 
a nice description of the process). The result is an inner 
transport of material and an outward transport of angular 
momentum. Numerical calculations of rotating magnetic 
disks indicate that a magnetorotational instability can ef- 
fectively yield a viscosity with a between 10~ 3 and 0.1 
(Brandenburg et al. 1999; Stone et al. 2000; Papaloizou & 
Nelson 2003). The suppression of the instability in neutral 
regions of the disk could imply rapid spatial variations of 
a (Gammie, 1996). 

Another mechanism initially proposed to be responsi- 
ble for angular momentum transport is thermal convection 
(Lin and Papaloizou, 1980; Ruden and Lin, 1986). The 
mechanism had largely fallen out of favor after direct hy- 
drodynamical simulations of thermal convection in disks 
by Stone and Balbus (1996; see also Stone et al. 2000) 
had resulted in very small angular momentum transport, 
mostly inward instead of outward. However, these results 
are challenged by recent simulations of convection in a 
baroclinically unstable disk that yield an outward angu- 
lar momentum transport and corresponding values of a 
between 10~ 4 and 10~ 2 (Klahr et al, 1999; Klahr and 
Bodcnheimcr, 2001). Turbulent angular momentum trans- 
port due solely to thermal convection can be tested by our 
models because of the particularity that it will cease in an 
optically thin medium. As proposed by Ruden and Pollack 
(1991), turbulent viscosity can then be modeled by a value 
of a which is uniform in the optically thick part, and goes 
to zero when the disk becomes optically thin. 

Another potential source of angular momentum trans- 
port is that due to shear instabilities (e.g. Dubrulle 1993). 
This mechanism was claimed to be inadequate on the ba- 
sis of numerical simulations (Balbus et al. 1996; Hawley 
et al. 1999). However, the disappearance of turbulence in 
the hydrodynamical simulations may be an artifact due 
to the limited resolution of the simulations (Longaretti, 
2002). Indeed, turbulence has been shown to be sustained 
in Couette- Taylor experiments with outward decreasing 
angular velocity (Richard and Zahn, 1999). A different 



prescription for the turbulent velocity may then be ap- 
plied: 



(4) 



where f2 is the disk rotation rate. Richard & Zahn found 
(3 ~ 2 x 10~ 5 . Hure et al. (2001) further showed that disk 
models using this prescription and (3 ~ 10~ 5 are as ca- 
pable of explaining the observed accretion rates and disks 
lifetimes as a models. A convenient feature of this turbu- 
lent parameterization is that v does not depend on the 
temperature: the disk's evolution then becomes indepen- 
dent of complex issues related to radiative transfer and 
opacities. 

Other mechanisms can yield angular momentum trans- 
port in circumstellar disks without necessarily being 
amendable to the calculation of a turbulent viscosity. This 
is for example the case of gravitational instabilities and 
density waves (e.g. Laughlin and Rozyckza, 1996; Laughlin 
et al. 1998), which may be an important source of disk 
evolution early on, when the disk/star ratio is still rela- 
tively large. It is also the case of bipolar outflows (e.g. 
Konigl, 1989; Casse and Ferreira, 2000). These proba- 
bly also have an important role during the cloud collapse 
phase. Although the ejection is limited to <~ 10% of the 
material accreting onto the central star (Calvet, 1997), 
this could nevertheless represent a substantial sink of an- 
gular momentum so that Eq. (2) would not be valid in the 
central region affected by the outflow. Since both of these 
effects are of importance only during the early evolution 
phases, neglecting them yields very limited quantitative 
changes to our results (see §5.1). 

3.3. Gravitational collapse of an isothermal molecular 
cloud core and disk growth 

The formation and early evolution of the star+disk system 
is governed by the collapse of its parent molecular cloud 
core. This process remains poorly known (e.g. Andre et al. 
2000), and we choose to model it using several simplifying 
assumptions. First, we assume that the cloud envelope is 
isothermal and spherically symmetric. Any given shell of 
radius I and angular velocity u(£) will collapse onto the 
disk within the centrifugal radius (the point at which the 
maximal angular momentum in the shell is equal to the 
angular momentum in the disk): 



R c (t) = 



i(t)Mi) 2 

GM(t) 



(5) 



where M(t) is the mass that has been accreted onto the 
star+disk system at time t. (In the formalism of the self- 
similar solutions, t = corresponds to the formation of 
the central condensation). Formally, eq. 5 is valid only in 
the limit when the disk's gravitational potential can be 
neglected. However, given our lack of knowledge of the 
collapse phase, this simplification is perfectly justified. 

We further assume that angular momentum is con- 
served during the collapse and that material falling onto 
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the disk finds its way to a location where its angular mo- 
mentum is equal to that of the (keplerian) disk. With that 
hypothesis, and given M, the accretion rate onto the disk, 
we derive the following source term: 



S(r,t) = 



M 1 



llR? 8 



-3/2 



1/2 



"I-V2 



(6) 



Departures from a central gravitational potential during 
the collapse phase imply that Eq. (2) does not perfectly 
conserve the angular momentum of the disk+envelope sys- 
tem. However, once again, the effect is negligeable. 

This expression of S(r, t) differs from that obtained 
from ballistic integrations by Cassen & Moosman (1981) 
who showed that the envelope material that encounters 
the disk has in fact a subkeplerian rotation rate (see also 
Nakamoto & Nakagawa 1994). This yields a small ad- 
ditional, outward, angular momentum transport that we 
choose to ignore in this simulation: as shown in Fig. 1 be- 
low, the effect would tend to decrease the surface density 
in the outer regions by 40% at most. The turbulent vis- 
cosities that are considered have a much more significant 
effect. Furthermore, this is also relatively small compared 
to most sources of uncertainties, in particular those related 
to the cloud collapse itself. 

The values of M and R c (t) depend on the mechanisms 
that led to the collapse of the molecular cloud. A simple 
and widely used solution is that of Shu (1977). In that 
case, a self-similar approach to the problem with uj(£) — 
Lo c d, a constant rotational speed of the molecular cloud, 
shows that the collapse may proceed from inside-out, with 
a constant mass accretion: 



M = 0.975-^, 



(7) 



where c s is the isothermal sound speed. Because the col- 
lapse solution yields £ = c s t/2 and using a mean molecular 
weight /i = 2.2, one can show that the centrifugal radius 
can be written: 

3 



y ' Vio-^s- 1 / 



Tod V 
10K J 



(M(t)Y 



AU, 



(8) 



Observations of hot cores indicate that the accretion 
rate is not constant over time but enhanced in the very 
first stages of cloud collapse and progressively diminish- 
ing after a more or less long time span of accretion rate 
close to the one predicted by the Shu theory (Bontemps et 
al. 1996). However, observations of some very young pro- 
tostars like IRAM 04191 in the Taurus molecular cloud 
(Belloche 2002) point to the presence of differential rota- 
tion and hence different initial conditions. Indeed, in the 
presence of a substantial magnetic field, the collapse phase 
can begin in a cloud whose rotation rate is u(£) oc 1/t. 
In that case, accretion proceeds much faster, and a self- 
similar approach shows that the centrifugal radius then 
evolves linearly with the accreted mass (Basu 1998): 



M~ 10 



G ' 



(9) 
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Fig. 1. Surface density obtained at the end of the collapse of 
a 0.3 Mo disk, assuming a centrifugal radius R c = 100 AU, 
and no viscous stress. The plain line corresponds to the col- 
lapse of a molecular cloud core in solid rotation. The dotted 
line shows the same solution when accounting for the outward 
angular momentum transport due to the subkeplerian momen- 
tum of the incoming envelope material (see Cassen & Moosman 
1981). The dashed line corresponds to the collapse of magne- 
tized cloud core with a rotation rate u)(£) oc l/£. 



R c (t) ~ 15 



( Y ( Bref \ 

VlO-^s- 1 / UomG/ 



(^) AU, 



(10) 



where cub is the ambient rotation rate of the cloud and 
B TC { is a background reference magnetic field. The numer- 
ical factor in front of the accretion rate is one possible 
value among many but all models imply it is of order 10 
(e.g. Foster & Chevalier 1993, Basu 1998, Hennebelle et 
al. 2003). Because the accretion rate also depends on the 
poorly- known c^, this uncertainty is ignored. 

In the absence of a viscous stress, it can be shown 
that the collapse of the cloud core yields a profile density 
E(r) = J S(r,t)dt that is approximately proportional to 
r~ 7 / 4 in the case of a cloud core in solid rotation and to 
3/2 m the case of our magnetized, differentially rotating 
cloud core. The resulting profiles are compared in Fig. 1. 
Because the centrifugal disk radius at the end of the col- 
lapse is fixed in this comparison, the magnetized cloud 
core has initially more angular momentum: the disk that 
forms therefore has more mass in its outer regions. 

3.4. Calculating disk properties 

Solving Eq. 2 necessitates the calculation of quantities 
such as the mid-plane temperature, T c , the vertical scale 
height, H, and the keplerian rotational frequency f2^. This 
has been carried out in many papers (see e.g. Ruden 
and Lin, 1986; Ruden and Pollack, 1991; Reyes-Ruiz and 
Stepinski, 1995 to cite only a few), but the detailed ex- 
pressions may differ slightly from one work to another. 
Generally, these works have considered disks of small 
masses and sizes, and have neglected the gravitational po- 
tential of the disks themselves. Because this effect may 
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be important for the evolution of DM Tau and GM Aur, 
we rederive here vertically averaged equations including 
autogravitation. 

3.4.1. Autogravitation and general disk properties 

We follow Ruden and Pollack (1991) in expressing the 
isothermal sound speed, c s , and the mid-plane density, p c 
as 



~ z — T 



E 

2H' 



(11) 



(12) 



Here 1Z is the gas constant, fi the molecular weight of the 
gas, supposed to be 2.2, and H is a measure of the vertical 
scale height, which is a well defined parameter only in non 
autogravitating disks. The mass of the disk also affects the 
keplerian rotational frequencies of the disk: 



f2 fe (r,i) = 



GM»(f) | ldV d 
r 3 r dr 



1/2 



(13) 



where Vd is the gravitational potential of the disk. This 
potential can be calculated by 



Vd(r) 



2tt 



-GS(n)n 



r, Jo \J (ri sin 9) 2 + (r — r[ cos 9) 2 



dnd6. (14) 



It is convenient to express Eq. (14) in terms of the 
radial coordinate alone by using the elliptic integral of 
the first kind K\m\: 



roc 4K 
V d (r) = / -G— 

JR. 



4(l+r/n-l)) 



\(r/n-l)\ 



■S(n)dn. 



(15) 



Generally, the incorporation of Vd does not play a ma- 
jor role in the evolution of the disk. Its effect is to slightly 
increase fife decreasing the turbulent viscosities through 
the definitions of H and v. Its influence is greater in mas- 
sive and extended disks with an inner low mass star. 

The vertical scale height, H, is calculated assuming 
vertical hydrostatic equilibrium, considering the vertical 
component of star gravity and the local gravitation of 
the disk. This is approximated from the infinite and R- 
homogeneous slab approximation, which considers the lo- 
cally equivalent limits S=cte or z — > (Paczyhski, 1978; 
Hure, 2000): 



ldP 

p dz 



= -Vl\z - AnGYi = gz, 



(16) 



where g is the effective gravity from both star and disk. 

An important simplification is to approximate the 
disk's vertical structure by an isotherm (the departures 
from an isotherm yield only second order effect). In that 
case, the disk pressure and density decrease exponentially: 



(P, p) = (Po, P o) cxp (-(z/Ho + (z/Hi) 2 )) 



(17) 



where the Hq term comes from the autogravitation term 
and Hi is the classical scale height coming from the verti- 
cal component of the radial gravity of the star. They are 
defined as 



Ho 
Hi 



4ttGE' 
V2C, 



(18) 
(19) 



The first one is dominant for large distances and large 
values of S, where the disk gravitation dominates over the 
vertical component of the central star gravitation. The sec- 
ond one is the scale height in absence of autogravitation. 
Consistency with our definition of p c through Eq. (12) 
requires that H has the following form: 



H = Hi — exp 



Hi 
2H 



1-Erf 



Hi 
2H 



(20) 



When autogravitation is negligible, Hq is large and H — > 
v / 7r/2ifi. When the disk's gravity is significant, Ho is 
small and H < Hi . The role of the disk gravity is mainly 
restricted to produce a slight vertical flattening of the 
outer disk. 

The assumption of vertical isothermal structure is gen- 
erally accurate only close to the disk mid-plane and values 
of H should be considered only as reasonable estimates 
instead of exact calculations. This factor yields a slight 
uncertainty on the magnitude of v in the case of the a 
parameterization. 

3.4.2. Temperature calculation 

In order to calculate mid-plane temperatures the following 
assumptions are made: 

— The disk is geometrically thin; 

— It is heated by the star's illumination but also by the 
dissipation of viscous energy by turbulence; 

— It is assumed to be optically thick in the radial di- 
rection everywhere so that heat can be transported 
efficiently only in the vertical direction. 

In thermal equilibrium, the disk's cooling flux is equal 
to the sum of the heating fluxes due to viscous dissipation 
and external sources: We define 7] as the effective temper- 
ature at which the disk is being heated by the central star 
and external sources and T e as the emission temperature 
at which the disk cools down. These quantities are then 
related by the following expression: 



2o- B T e 4 = Y,vr 2 tt 2 r + 2a B T l 4 , 



(21) 



where <jb is Stephan-Boltzmann's constant, T e is the effec- 
tive temperature of the disk, and 7] is an effective temper- 
ature to which an inert disk would be heated by external 
sources. T; will be discussed in the next section. 

Obtaining the mid-plane temperature T m from the ef- 
fective temperature T e is generally a complex radiative 
transfer problem. However, in both the optically-thick 
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and optically-thin regimes, analytic expressions can be 
found (assuming that the opacities behave correctly). In 
the present work, we adopt the expression derived by 
Nakamoto and Nakagawa (1994), who approximate the 
mid-plane temperature as a sum of optically-thick and 
optically-thin contributions: 



T4 1 ( 3 1 



Y,vr 2 nl + a B T, 



i > 



(22) 



where tr and Tp are the Rosseland and Planck mean op- 
tical depths respectively. The optical depth from the mid- 
plane to the surface of the nebula, tr, is defined in terms 
of the Rosseland mean opacity K R (p c ,T c ) by 



tr = 



(23) 



To evaluate kr, we used the Rosseland mean opacity of 
dust grains provided by Pollack et al (1986) in the form of 
a power law: kr(p c ,T c ) = K p c m T c n . Note that the rapid 
changes of the opacities in different regimes (e.g. due to 
the evaporation of water, ...etc.) imply that a robust au- 
toconsistent numerical method must be used to solve the 
temperature equation. More recent opacities provided by 
Hartmann and Kenyon (1996) and Bell and Lin (1994) 
seem to be in general agreement with the values of Pollack 
et al. As Nakamoto and Nakagawa (1994), we further as- 
sumed Tp = IAtr. 

Since the intensity of the turbulence viscosity is de- 
termined by the disk temperature and determines itself 
part of the heating the set of nonlinear equations must be 
solved simultaneously in an autoconsistent manner once 
a value of a is given. A time-forward integration of the 
density surface distribution is then possible. The problem 
is more straightforward in the case of the (3 prescription 
because v is then independent of temperature. 

3.4.3. Irradiation from the central star 

Stellar photons reprocessed by the accretion disk consti- 
tute a significant contribution to the disk's thermal bud- 
get, especially in regions where dissipation by viscous pro- 
cesses is small, i.e. mostly in the external regions of the 
disk. Again, a proper treatment of stellar irradiation is 
beyond the scope of this work, but we attempt to capture 
the essential underlying physics under our ID approach. 

At distances much larger than the stellar radius, it 
can be shown that stellar irradiation implies an effective 
temperature that is a function of the star's effective tem- 
perature T* and radius i?* (Adams et al., 1988; Ruden 
and Pollack, 1991): 



3tt V r ) 



1_,R. 
2 V r 



) 3 (f)(^f-0 



1/4 



(24) 



The first term inside the brackets is the contribution due 
to the finite-sized star irradiating a flat disk. The second 
term, proportional to (d In H/d In r — 1) accounts for the 
flaring of the disk. We find numerically that this factor 



becomes significant only for radii r > 50 AU. However, 
this factor is a significant source of numerical instabili- 
ties. This is an interesting problem that certainly requires 
further work. At present, we choose to avoid it by impos- 
ing 

d\nH/d\nr = 9/7, 

which corresponds to the (approximate) equilibrium so- 
lution for a disk whose temperature is dominated by the 
flaring term (Chiang & Goldreich 1997). We verified that 
this hypothesis is autoconsistent, i.e. our disks are close 
to H oc r 9 / 7 in regions where the flaring term dominates. 

We additionally considered that the disk is heated by 
its environment which radiates at the same temperature 
as the molecular cloud, T; 2 = T C( j. The true irradiation 
temperature T\ is then computed as 

T?=T l 4 i +T* (25) 

3.5. Comparison to an a-disk model with detailed 
radiative transfer. 

Our model was tested against a calculation by d'Alessio 
et al (1999, 2001) for a 0.5 M Q T-Tauri star. Their calcu- 
lation of the density profile is based on a standard a-disk 
model, in which a static solution to Eq. (2) with no source 
term is found by imposing a constant mass flux throughout 
the disk. As shown in fig. 2, the resulting surface densities 
for a given accretion rate onto the star differ slightly. This 
is not unexpected, due to the fact that we included the 
collapse phase, because assumed opacities are different, 
and because the static solution is not necessarily a good 
estimate of the disk's structure, especially at large orbital 
distances. 

Despite these differences, the midplane temperatures 
inferred in our model compare well to the more elaborate 
2D calculations by D'Alessio et al. Figure 3 shows that the 
comparison is good in two regimes: (i) In the inner parts of 
the disk (R 10 AU), where the release of gravitational 
energy by viscous dissipation dominates over the star's 
irradiation and governs the disk's temperature structure; 
(ii) In the intermediate and outer parts of the disk, where 
irradiation from the central star is the dominant source of 
heating and the temperature profile tends towards T(r) oc 

r -l/2_ 

3.6. Limitations of the model 
3.6.1. Structure of the inner disk 

In the calculations that are presented, we assume that the 
disk extends from an arbitrary inner boundary at 0.1 AU 
and beyond. The exact location of this inner boundary is 
of little importance for the evolution of the outer disk (say, 
beyond 10 AU). However, the exact structure of the inner 
disk and the fact that recent observations seem to point 
to a clearing of the region inside <~ 5 AU (Bergin et al. 
2004) can affect the constraints derived from the accretion 
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0.1 1 10 100 

R [AU] 

Fig. 2. Surface density as a function of orbital distance for 
three different times, corresponding to accretion rates onto the 
central protostar of 10~ 7 , 10~~ 8 and 10 -9 Mq yr -1 . Plain lines 
show the results of our model. These are compared to similar 
calculations by D'Alessio et al. (2001) for a standard a-disk 
with a = 0.01, M, = 0.5 M , R, = 0.0093 AU and T, = 
4000 K (dotted lines). We further used: M = 0.2 M Q , ui cd = 
3 x 10~ 3 s _1 and T c d = 10 K implying a total accretion time of 
the cloud core of 0.3 Myr and an outer maximum centrifugal 
radius of 60 AU. 




10 r 



0.1 1 10 100 

R [AU] 

Fig. 3. Midplane temperatures as a function of orbital distance 
calculated in this work (plain lines) and by D'Alessio et al. 
(2001) (dotted lines), in the same conditions as in Fig. 2. 



rate. For example, it can be noticed that the presence of 
a Jupiter-mass planet would lower the accretion rate onto 
the star, because part of the material from the outer disk 
would be accreted by the planet (e.g. D'Angelo et al. 2002, 
Bate et al. 2003), and part of the disk's gravitational en- 
ergy would be used to push the planet inward. Clearly, in 
the absence of detailed models of this and because of other 
possible explanations (photocvaporation, grain growth) it 
would be unrealistic to try to also fit the observed struc- 
ture of the inner disk with our simple model. Last but 
not least, the error bar on M is relatively large and is 
presumably larger than this source of uncertainty. 



3.6.2. Thermal vertical structure 

In models with an a- viscosity, the temperature determines 
the sound speed c s , the scale height H and hence the vis- 
cosity v itself. We use the midplane temperature because 
the gravity near the midplane is low (~ uz) and vertical 
temperature variations should be small in these regions. 
However, radiative transfer models (Malbet & Bertout 
1991; Chiang & Goldreich 1997; D'Alessio et al, 1999, 
2001) show that in regions where the stellar irradiation 
is significant, the temperature in the outer layers of the 
disk is expected to be larger than at the midplane. This 
was verified by (Dartois et al., 2003) who found that in 
DM Tau's outer disk, the central temperature is T c ~ 15 K 
while Tphotosphoro ~ 30 K. 

Another aspect of the problem is that in the tenu- 
ous atmosphere of the outer disk, the gas may not be in 
thermal equilibrium and may therefore be heated more 
than the dust disk, i.e. by absorption of FUV photons. 
Altogether, this should introduce a maximum of a factor 
of 2 indetermination on the temperatures and a- viscosities 
of the outer optically thin disk. 

3.6.3. Non-viscous processes affecting the disk 
evolution 

Our model assumes that the present disks around DM 
Tau and GM Aur are solely the result of the collapse of 
a molecular cloud core and of the subsequent viscous evo- 
lution of the disk. May these disks be altered by other 
physical processes? Following Hollenbach et al. (2000), we 
identify five possible sources of disk dispersal: 

— Wind stripping; 

— Photoevaporation due to an external source; 

— Photoevaporation due to the central star; 

— Tidal stripping due to close stellar encounters; 

— Planet formation. 

Wind stripping may indeed affect the disk thickness, 
but according to Hollenbach et al. (2000), the associated 
timescales for disk dispersal are, in the regions of interest, 
well over 10 Myr and can be neglected. 

Unlike in Orion (e.g. O'Dell et al., 1993), photoevapo- 
ration due to external sources should not be a concern in 
the Taurus- Auriga association due to the small number of 
massive stars present. 

Photoevaporation due to the central star may be esti- 
mated from the work of Hollenbach et al. (1994; see also 
Shu et al. 1993). The authors estimate that the total mass 
loss due to UV photoevaporation of disk material is: 

M uv = 4.1 x lO" 10 ^ 2 M* 1/2 M yr-\ (26) 

where $41 is the ionizing photon luminosity of the central 
star in units of 10 41 s _1 . There is no precise value of this 
parameter, but Alexander et al. (2005) advocate using the 
observed He II/C IV line ratio to derive the hardness of the 
ultraviolet spectrum and obtain $41 ml — 1000 for five 
classical T-Tauri stars. As shown by Clarke et al. (2001) 
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and Matsuyama et al. (2003), photoevaporation alters sig- 
nificantly the inner disk's structure when Muv ~ M* 7 
where M* is the accretion rate onto the star. This implies 
that photoevaporation may play a role even in DM Tau 
and GM Aur if the UV photon rate is in the upper range 
of that estimated by Alexander et al. Even in that case 
however, we expect the effect on the structure of the outer 
disks to be relatively modest. 

Tidal stripping may be important for the evolution of 
accretion disks in relatively dense environments. A simple 
estimate for the timescale of truncation of a disk to a ra- 
dius rd assumes that the disk is stripped to about 1/3 the 
impact parameter (Clarke and Pringle, 1991, Hollcnbach 
et al. 2000); hence: isE ~ l/(n*av), where is the den- 
sity of stars, a ~ ir(?>rd) 2 is the collision cross section and 
v is the velocity dispersion of stars. In our case, assuming 
tsE ^ 10 7 yr, v i=a lkms" 1 and r<j « 1000 AU implies that 
stellar encounters could become important for star den- 
sities n± > 50 pc~ 3 . Values n± w 100 pc~ 3 appear to be 
typical of the central Taurus cloud (Clarke et al. 2000). 
Moreover, both DM Tau and GM Aur are relatively iso- 
lated in this cloud and consequently this mechanism has 
probably not affected their disks. This is probably unlike 
other young stars, which may explain why young stars 
with very extended disks remain relatively rare. 

Last but not least, planet formation is a mechanism 
susceptible of greatly modifying the structure of disks, 
particularly when giant planets form. However, planets 
are thought to form closer to the central star than can be 
probed by the present observational techniques. Their ef- 
fect on the disks can be thought as a surface density sink, 
and should hence be relatively modest at large orbital dis- 
tances. 



3.6.4. Gravitational stability of the disk 

The local gravitational stability of a rotating disk against 
axisymmetric perturbations is measured by the Toomre-Q 
parameter, which is defined by 



Q 



kc s 



(27) 



where k is the epicyclic frequency given by k 2 = 
{l/r 3 )[d{r 4 fl k 2 )/dr} (Toomrc, 1964; Goldrcich and 
Lynden-Bell, 1965). Disks are gravitationally stable if 
Q > 1 all over the disk (the rotation gradient domi- 
nates and the disk follows axisymmetric evolution) and 
are unstable if Q < 1 anywhere (gravitation may domi- 
nate locally and break out the axisymmetry) . In this case 
disks may produce a direct gravitational collapse of parti- 
cles into planetesimals (Goldreich and Ward, 1973), giant 
planet cores (Boss, 2003) or may develop spiral density 
waves able to redistribute the angular momentum in a 
non local way until stability is reached again (Laughlin 
and Rozyczka, 1996; Laughlin et al. 1997, 1998). This last 
case is the generally accepted scenario for the evolution of 
massive disks. 



Nakamoto and Nakagawa (1994) showed that disk in- 
stabilities are more likely to appear in disks with low val- 
ues of a and high values of kJ c d- Laughlin and Rozyczka 
(1996) found that in most of the situations the global 
transport of angular momentum by spiral arms in gravita- 
tionally unstable disks may be roughly in agreement with 
axisymmetric models with values of a s of 0.02 — 0.03. 

Although both DM Tau and GM Aur disks are at 
present gravitationally stable (Q > 1 everywhere), they 
may have gone through an early unstable phase. In order 
to account for this effect, we adopted a = Max(0.03, a) 
when the condition Q < 1 was met somewhere in the disk. 
A limited number of models with a > 0.03 went through 
an unstable phase without any artificial increase in their 
viscosity. We did not modify the viscosity of /3-models, 
whatever the value of Q. 

We observed a posteriori that the gravitational insta- 
bility phase was generally limited in time to less than 
0.1 Myr. An incorrect handling of gravitational instabili- 
ties hence yields a negligible error on the age of the system. 
An effect that is potentially more significant, and ignored 
in the present study, concerns the fact that the advec- 
tion of disk material can depart from a diffusive solution. 
To the extreme, this process may lead to the formation of 
stellar or planetary companions. Clearly, a more consistent 
treatment of gravitational instabilities would be desirable. 

4. Fitting the observations 

4.1. Setting of the calculations 

The equations described in the previous section are solved 
numerically using an explicit finite-difference scheme. 
Each calculation begins at time t — Mq/M. The inner 
disk boundary is set to Ri n =0.1 AU, and the maximum 
computational space ends at R ou t — 10 4 AU. In order to 
solve numerically Eq. (2) we follow the method described 
by Bath and Pringle (1981). This requires the radial points 
to be equally spaced on a X = -Jr space. The stability 
condition required to solve Eq. (2) is that the time-step 
At obeys 



At > Min 



/X 2 AX 2 



24j/ 



(28) 



An adaptive time-step scheme taking into account this 
condition is used. A grid resolution of 250 points is ade- 
quate for disks with a centrifugal radius that is always 
larger than 4AU. To calculate the evolution of disks with 
smaller values of j c d so that 2 < R c < 4 AU, we use 500 
points, corresponding to a resolution in the inner regions 
- 0.15 UA. 

The model parameters are: 

— a or j3, which characterize the magnitude and func- 
tional form of the turbulent viscosity; 

— cj c d, the angular velocity of the molecular cloud core; 

— T c d, the temperature of the molecular cloud (which is 
equivalent to a characteristic accretion velocity from 
the molecular cloud onto the circumstellar disk); 
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Table 2. Explored parameter space 





DM Tau 


GM Aur 


Mo [M 1 


0.05 and 0.30 


0.05 and 0.40 


M cd [M ] 
tj c d [s _1 ] 


0.4 - 1.0 


0.4 - 1.5 


10~ 15 - 10~ n 


io- 15 - io- 11 


T cd [K] 


3.0-30.0 


3.0 - 30.0 


Accretion Timc(*) [yr] 


45000-5 xlO 6 


49000- 5 xlO 6 


a 


10~ 5 - 1 


10~ 5 - 1 


P 


10~ 6 -0.1 


10~ 6 -0.1 


Number of models 


2000 a 


2000 a 




1500 P 


1500 (3 


Magnetic cloud core 


2000 a 





(*) Here the accretion time refers to the time it takes for the 
molecular cloud core to entirely collapse onto the disk (see 
Eq. 7). 



— Mo, the initial mass of the protostar; 

— M c d, the total amount of material initially in the 
molecular cloud and eventually in the star+disk sys- 
tem. 

These parameters reflect either uncertainties in the 
physics of star formation (a- (3, Mo) or unknown initial 
conditions (w c( j, M c( j), or both (T c( j through the accretion 
rate of the molecular cloud, see Eq. 7) . 

Our approach is to constrain these parameters by ex- 
tensive numerical simulations and comparison with obser- 
vations. Four sets of calculations are performed for DM 
Tau and GM Aur, using either the a or (3 parameteriza- 
tions and considering the Shu (1977) classical cloud col- 
lapse model. To explore the sensitivity of the results to 
the cloud collapse model, we calculate a fifth set of model 
for the DM Tau a case and the collapse scenario derived 
from Basu (1998). 

In order to efficiently explore the space of parameters 
we use a Monte-Carlo procedure to chose the values of 
the model parameters within the range given in table 2. 
Two values of Mo are chosen. A small value is representa- 
tive of cases in which the central condensation grows rel- 
atively slowly, and gravitationnal instabilities dominate 
the early disk evolution. A large value is representative 
of an initially evolved situation in which the proto-star 
forms rapidly, and the disk grows significantly only at later 
times. This could occur if ambipolar diffusion and jets al- 
low an efficient outward transport of angular momentum 
early on during the collapse of the molecular cloud. 

Cases for which the initial centrifugal radius R c (see 
Eq. (5)) is larger than ^5000 AU are disregarded, as well 
as those with final centrifugal radius too small to solve 
accurately the building-up of the disk. For numerical rea- 
sons models that require a time-integration with a step 
smaller than 0.5 yrs are not calculated. Table 2 provides 
a summary of the different sets of calculations run and 
the coverage of the space of parameters. Table 3 pro- 
vides the numerical limitations imposed by the two dif- 
ferent grid resolutions, the minimum time-step required, 



Table 3. Numerical limitations of the calculation 



Npoints — 250 


Npoints =500 


AR = 0.4 


ARo = 0.15 


2.3 < R c < 5000 


0.8 < R c < 5000 


19.0 < Logio 0'cd) < 21 


18.8 < Logio (jed) < 21 


At > 0.5 => 




Logio M < 17.2 


Logio (vo) < 16.5 



ARo is the spatial resolution of the model at the inner disk 
boundary and R c is the centrifugal radius, both are given in 
AU, jed is the specific angular momentum of the molecular 
cloud and is given in cm 2 s _1 , At is the minimum allowed 
computational time-step in yr and is related with vo the value 
of the disk viscosity at the inner boundary given in cm -2 s _1 . 
The inner disk boundary is located at 0.1 AU. 

and the resulting limits on j c( j and i/q- The values of j c d 
bracket inferred values for Class I protostellar envelopes, 
10 20 < jed ^ lO^cm^- 1 (Ohashi et al. 1997). The up- 
per limit on v a is relatively large and limits us only for 
extreme values of the /3-parameter in the case of DM Tau 
(see figure 14 hereafter). 

The calculations were performed on an 8PCs 
Linux/Beowulf cluster over several months. The low At 
required by the stability criterion imply total computation 
times for each model that range between minutes to sev- 
eral hours for the most computationally expensive models 
(specifically, these correspond to 500-points, (3 models of 
GM Aur). The large database of resulting models is then 
compared to the various observational constraints in order 
to extract the possible values of the main physical param- 
eters. 

4.2. Selecting models from observational constraints 

The observational constraints characterizing DM Tau and 
GM Aur have very different origins. Three independent 
constraints are used (see § 2): 

(a) The fact that the disks are optically thick in 12 CO at 
100 AU and the outer disk radius requires that S be 
larger than a minimal value. This value is calculated 
by assuming that all C is in form of vapor CO and 
translates into a minimal value of £ at R ou t ■ 

(b) Values of £ at 100 AU have to be consistent with those 
inferred from dust emission, allowing for uncertainties 
in the opacity coefficient, grain size and dust to gas 
ratio. 

(c) The accretion rate onto the star has to fit constraints 
obtained from the star's excess luminosity. A conser- 
vative error bar of 0.54 in Logio (M) is assumed. 

In order to further narrow the possible solution ensemble, 
we choose to add the following reasonable constraints: 

(d) At the outer 12 CO radius, the value of S cannot be 
larger than 20 times the lower limit. This value is in- 



12 



Hueso and Guillot: Protoplanetary disks evolution: DM Tau and GM Aur 




® 

1 *■ logR 

100 AU outer disk 

Fig. 4. Illustration showing the different observational con- 
straints used in this work, in terms of surface density and or- 
bital distance. The labels "c" and "e" correspond to constraints 
on the accretion rate onto the central star. Constraints "e" and 
"d" are reasonable expectations of the mass accretion rate and 
outer disk maximum surface density. 

spired by the discrepancy between the CO and dust 
emission observations at 100 AU. 
(e) The inferred accretion rate is supposed to be known 
within a realistic but model dependent error bar of 
only 0.3 in Logio(M). 

These constraints are schematically shown on Figure 4. 
They should be further refined by combining observational 
and theoretical studies. In this work, we purposely use 
rather unrestrictive interpretations of the observations to 
obtain robust constraints on the desired physical param- 
eters. 

Together with the constraints on the star's mass and 
age, we thus derive five sets of solutions, from weaker to 
stronger constraints: 



{1} 


(a) 








{2} 


(a)- 


-(b) 






{3} 


(a)- 


-(b). 


-(c) 




{4} 


(a)- 


h(b)H 


-(c)- 


Kd) 


{5} 


(a)- 


-(b)- 


-(c)- 


h(d)+(e) 



We will mostly discuss the results obtained with sets {4} 
and {5}. 

4.3. Examples of models behavior: Disk history and 
surface density-temperature evolution 

Before analyzing our global results, it is useful to describe 
the general behavior of the models in terms of two spe- 
cific examples. We explain here how the observational con- 
straints are used to select a model as a successful fit to the 
data. We consider two examples based on an a-model of 
DM Tau. The parameters for both examples are listed in 
Table 4 

In all cases, our story starts with a protostar of mass 
Mo, no disk and a reservoir of cloud material with mass 
M c( j— Mo. Material falls from the molecular cloud at a con- 
stant rate inside a disk of growing size R c , a consequence 



Table 4. Model parameters for examples 1 and 2 



Parameters 

(fixed) 


Example 1 


Example 2 


Units 


a 


0.01 


0.025 






2.3 x 10" 14 


2.6 x 10" 13 




T c d 


14 


17 


(K) 


M 


0.05 


0.05 


(M ) 


M cd 


0.515 


0.585 


(M s ) 


(derived) 


Logio(Jcd) 


52.3 


53.4 


(gcm 2 s 1 


Logio(jcd) 


19.3 


20.3 


(em's- 1 ) 


R c 


11 


830 


(AU) 


Accretion time 


0.18 


0.15 


(Myr) 



Star Age 
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/ /' 
/ / 


''Mac rc 
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Fig. 5. Example 1. Evolution of star mass M* and disk mass 
Afdisk as a function of time with masses in solar units (corre- 
sponding axis to the left) for the model parameters of Table 4. 
The accretion rate onto the central star is shown as a dot- 
ted line (corresponding axis: first to the right). The evolution 
of the centrifugal radius R z (see Eq. (5)) is shown as a dash- 
dotted line (corresponding axis: far-right). Gray curves and the 
hashed region indicate time sequences when selected observa- 
tional constraints are verified (see text). 

of angular momentum conservation and the inside-out col- 
lapse of a molecular cloud core in solid rotation. Because 
viscous diffusion is more efficient than the increase of R c 
with time, the disk spreads beyond the centrifugal radius. 
The disk can be quite hot in this early phase, especially 
when R c is small and the accretion rate is large (large 
T cc j). The disk grows in mass until all of the available 
mass in the molecular cloud core has been accreted to the 
disk. After that, the central star continues to accrete disk 
material while the disk expands and cools. 

Starting with Example 1, fig. 5 illustrates the evolu- 
tion of the star mass, disk mass, centrifugal radius and 
star accretion rate. Figure 6 shows the evolution of the 
surface density of the disk. The disk construction phase 
(dotted lines) appears as a quick phase in which the disk 
density increases by orders of magnitude with a relatively 
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Fig. 6. Example 1. Surface density versus orbital distance at 
different times for the model parameters of Table 4. Dotted 
lines correspond to the early formation of the disk. The collapse 
of the molecular cloud ends after 0.18 Myr. The dashed line at 
10 Myr corresponds to the end of the simulation. The error bars 
represent the E values at 100 AU and in the outer radius that 
are used as observational constraints for DM Tau. The gray 
area shows the ensemble of models fitting those constraints. 
The dark-shaded region shows the E distribution at the time- 
lapse when all observational constraints are satisfied (see text). 




10 I I 

0.1 1 10 100 1000 10000 

R [AU] 

Fig. 7. Example 1. Same as fig. 6, but for the midplane tem- 
perature as a function of orbital distance. 

sharp outer edge. When the molecular cloud has collapsed 
entirely (at t = 0.18 Myr), the disk has already spread 
further than the maximum centrifugal radius (11 AU), a 
consequence of the relatively large v. The later evolution 
of the system is characterized by the inner accretion of 
disk material and the outer diffusion of angular momen- 
tum. After 3 Myr, the diffusion of material proceeds much 
more slowly due to the lower density, temperature and vis- 
cosity of the inner disk. 

The midplane temperatures are shown in Fig. 7 for 
the same time-sequence. The disk heats up considerably 
during the molecular cloud collapse phase (the first 0.2 
Myrs) and cools down gradually when the collapse ceases. 
The sharp transitions in the temperature curves arise from 
the different regimes of Rosseland opacities at different 
temperatures (Ruden and Pollack, 1991). The outer disk 
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Fig. 8. Example 2. Evolution of star mass M* and disk mass 
M disk as a function of time. See fig. 5 and Table 4 for details. 

beyond 50 AU is optically thin, vertically extended, flared, 
and its thermal structure is determined solely by stellar 
irradiation and a diffuse heat source of T cc i- 

In Example 1, the disk is always stable to gravitational 
perturbations (Q > 1) and the disk mass is never a large 
fraction of the star mass. This model evolves smoothly 
with values of the initial parameters in good agreement 
with expected values in the Taurus Aurigae region. The 
question then is: Does this model satisfy the observational 
requirements for DM Tau? And in that case, what other 
values of the set of parameters representative perhaps of 
very different initial conditions or turbulence in the disk, 
would also agree with the observations? 

Figure 5 shows thick grey lines superimposed on each 
plotted quantity. Each one shows the range of time for 
which that quantity agrees with the available observations. 
For Moisk, the grey line represents the period of time 
when the £ surface density satisfies the observational error 
bars discussed in Section 2. This period is easily identified 
on Figure 6. The accretion rate is reproduced either within 
large error bars (thick grey line), or with small error bars 
(thick black line) (see Section 2). The uncertainty over 
the star age for DM Tau is marked as a light-grey box. 
Figure 5 hence shows that the model is a good fit to the 
data from 1.5 to 2.8 Myrs (surface densities and star age), 
from 1.5 to 2.6 Myrs [M with its large error) or from 1.5 
to 1.6 Myrs (small error bars on M). The latter is shown 
as a dashed region. 

This model hence does fulfill the "strict" observational 
constraints (set {3}) and also the reasonable additions 
(sets {4} and {5}). This example shows how DM Tau's 
800 AU disk can be formed by viscous diffusion of an ini- 
tially much smaller disk, with a centrifugal radius i? c =ll 
AU. 

Let us now examine our second example. Here the disk 
is assumed to form in 0.15 Myrs with a maximum centrifu- 
gal radius of 800 AU. Most of the disk material falls so far 
from the central star that the disk gets more massive than 
the central star at the end of the collapse (figs. 8 and 9). 
Yet, the disk diffuses outwards and gets accreted into the 
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Fig. 9. Example 2. Surface density versus orbital distance at 
different times. See fig. 6 and Table 4 for details. 
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Fig. 10. Example 2. Midplane temperature versus orbital dis- 
tance. See fig. 7 and Table 4 for details. 



central star. For a certain period of time (4-5.4 Myrs) it 
also satisfies all the observational constraints wc have con- 
sidered. 

Observationally it may be very difficult to distinguish 
between these two scenarios at a late and evolved age but 
each one has very different implications for the formation 
of planetary systems. In the first case, a massive, rela- 
tively small disk forms. This may allow a rapid build-up 
of planetesimals and protoplanetary cores, but may induce 
a fast inward migration. In the second example, planetes- 
imals could still be formed in an extended disk on longer 
timescales. If massive cores were to be formed late in the 
system their inner migration could be damped. The possi- 
bility of forming planetesimals and planets in this variety 
of scenarios will be studied in a forthcoming paper. 

As we will see in the next section, solutions matching 
the observed surface densities and other constraints are far 
from being limited to these two cases; a variety of initial 
conditions and values of the turbulent viscosity can yield 
satisfactory models. We now turn to an analysis of the 
ensemble of models that match the observations. 



5. Results: Allowed models and derived 
theoretical constraints 

5.1. General 

For the most of this section we discuss results on the basis 
of only 2 parameters representative of the disk physics (a 
or (3) and of the initial conditions {j c d), respectively. We 
thus derive constraints on these parameters depending on 
the chosen set of observational constraints. Since given 
values of (a or (3, j c( j) may be obtained from different 
combinations of the initial conditions, we also plot the 
fraction of models that fit the observations. We further 
checked that the global results are not dependent on the 
value of Mq or on the numerical resolution. 

Figure 11 constitutes a summary of our global re- 
sults for DM Tau and GM Aur for both parameteriza- 
tions of turbulence. Results are shown for three sets of 
observational constraints ({3}, {4} and {5}). The most 
extended light-grey contours correspond to the envelope 
in the (v,j c d) space covered by models matching the CO 
and dust observations, and the star accretion rate with the 
large error bar (set {3}). This results in relatively weak 
constraints on the physical parameters. Less restrictive ob- 
servational constraints are therefore not shown: For DM 
Tau a models, 30% of all launched models satisfy the CO 
data alone (set {1}) and 18% satisfy the CO+dust data 
(set {2}). If we include the accretion rate information we 
are able to slightly reduce the number of fitting models 
to 15% (set {3}). These constraints are quite unrestrictive 
because they do not limit the outer radius of the disk and 
unreasonably massive and extended disks can result. 

Adding a reasonable maximum value for the density of 
the disks at their outer edge (set {4}) yields a significant 
reduction of the permitted (v, j c( j) space, either for the DM 
Tau (3 model or for the GM Aur a and (3 models. This is 
shown as a dashed contour and grey color in fig. 11. A still 
narrower set of solutions can be achieved by restricting 
the range of allowed accretion rates (set {5}). The final 
ensemble of solutions is plotted as a thick colored contour 
in fig. 11. We consider that this set of solutions provides 
the most realistic constraints on v and j cc i for DM Tau 
and GM Aur. 

Figure 11 also shows that a right value of a, (3 and 
jcd is not a sufficient condition to match the observations. 
The presence of other initial parameters imply that only 
a fraction of the models really fit the imposed conditions. 
This fraction peaks to 60% for the central region in the 
DM Tau a panel, 40% for the j3 case and 50% and 30% 
for GM Aur a and (3 respectively. 

The detailed values of the space of parameters allowed 
for sets {2} to {5} are given on Table 5. We list the max- 
imum, minimum and statistic mean value of the model 
parameters that fit each set of conditions. 

Our first conclusion is that it is possible to use the 
estimated current state of observed disks to get informa- 
tion about their formation and evolution using relatively 
simple models. Globally the constraints on the parame- 
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Fig. 11. Viscosity and initial specific angular momentum of models fitting DM Tau and GM Aur for both parameterizations 
of turbulence. The contour colored plot represents the fraction of models fitting the strongest observational constraints (set 
{5}). The light-shaded region around, contoured by a dotted line, represents the additional area covered by models assuming 
an uncertain accretion rate (set {4}). The light grey filled region surrounding this area corresponds to models with no upper 
limit on the surface density in the disk's outer region (set {3}). The spatial resolution for each panel is shown as an empty circle. 
The circle in each figure shows the approximate spatial resolution of the figure obtained by the Monte-Carlo procedure for an 
accuracy in the fraction of models of 10%. Features in the plots smaller than each circle are not well resolved. 



ters are not very strong, but they anyway provide useful 
information and ask for further detailed observations of 
these objects. For instance, a reduction of the star age un- 
certainty and/or the star accretion rate error bar would 
be especially useful. Another important conclusion to be 
drawn is that both a and f3 models provide a satisfac- 
tory evolutionary scenario for the disks around DM Tau 
and GM Aur. This does not mean that either of these pa- 



rameterizations are correct but that none of them can be 
excluded on the basis of these results. 

The current disks around DM Tau and GM Aur can be 
obtained from very different values of the angular momen- 
tum in the molecular cloud, corresponding to values of the 
centrifugal radius spanning the entire range (4 to 5000 AU) 
allowed by our numerical approach. A very small value of 
R c requires large values of the viscosity to yield an efficient 
outward diffusion of the disk to match its outer radius. A 



16 



Hueso and Guillot: Protoplanetary disks evolution: DM Tau and GM Aur 



System Param. Model results obtained for the different observational constraints 







{2}: CO + Dust 


{3}: {2} + M 




{4}: {3} + Outer Limit 


{5}: {4} + Accurate M 






Min 


Mean 


Max 


Min 


Mean 


Max 


Min 


Mean 


Max 


Min 


Mean 


Max 


DM Tau 


a 


5xlCT 4 


0.03 


0.6 


5xl0~ 4 


0.02 


0.15 


5xl0~ 4 


0.02 


0.15 


8xl0~ 4 


0.02 


0.08 




jcd 


19.0 


20.1 


20.8 


19.2 


20.1 


20.7 


19.2 


20.0 


20.7 


19.3 


20.0 


20.7 




Rc 


2.4 


1090 


4900 


4.6 


980 


4900 


4.6 


720 


4900 


7.3 


750 


4900 




<^cd 


0.14 


22 


218 


0.14 


23 


218 


0.14 


22 


218 


0.14 


23 


218 




Ted 


2 


15 


30 


2 


16 


30 


2 


16 


30 


2 


16 


30 




M cd 


0.41 


0.59 


1.0 


0.41 


0.56 


0.87 


0.41 


0.53 


0.69 


0.41 


0.53 


0.69 




Jcd 


52.0 


53.5 


54.1 


52.2 


53.1 


54.0 


52.2 


53.0 


53.9 


52.3 


53.0 


53.7 




AT 


0.05 


0.5 


4.4 


0.05 


0.4 


4.0 


0.05 


0.3 


4.0 


0.05 


0.3 


3.7 




Age 


1.5 


3.7 


7.0 


1.5 


3.9 


7.0 


1.5 


3.5 


7.0 


1.5 


3.5 


6.8 




AAge 


0.05 


2.6 


5.5 


0.05 


2.3 


5.5 


0.05 


2.0 


5.5 


0.05 


1.0 


2.8 




M disfe 


0.006 


0.05 


0.42 


0.009 


0.036 


0.38 


0.009 


0.025 


0.10 


0.014 


0.036 


0.097 




^models 




371 






313 






225 






166 




DM Tau 


/3xlO b 


1.1 


31 


390 


1.1 


24 


160 


1.1 


8.5 


50 


2.2 


8.5 


50 




Jcd 


19.1 


20.0 


20.8 


19.3 


20.1 


20.8 


19.3 


19.7 


20.1 


19.4 


19.7 


20.0 




Rc 


5 


752 


4800 


9 


810 


4800 


9 


74 


328 


15 


71 


170 




^cd 


0.05 


14 


140 


0.1 


15 


140 


0.1 


8 


36 


0.1 


8 


36 




Ted 


2 


14 


30 


2 


15 


30 


2 


15 


30 


2 


15 


30 




M cd 


0.413 


0.65 


0.97 


0.42 


0.65 


0.97 


0.42 


0.56 


0.72 


0.43 


0.55 


0.65 




Jcd 


52.1 


53.1 


54.1 


52.4 


53.2 


54.1 


52.4 


52.8 


53.1 


52.5 


52.8 


53.1 




AT 


0.06 


0.7 


4.9 


0.06 


0.6 


4.9 


0.06 


0.4 


2.6 


0.06 


0.5 


2.6 




Age 


1.5 


4.0 


7.0 


1.5 


4.0 


7.0 


1.5 


3.2 


7.0 


1.5 


2.7 


4.5 




AAge 


0.05 


2.7 


5.5 


0.05 


2.0 


5.5 


0.05 


0.8 


3.4 


0.05 


0.4 


1.3 




Mdi s fc 


0.01 


0.15 


0.57 


0.01 


0.20 


0.50 


0.01 


0.06 


0.13 


0.02 


0.07 


0.13 




N models 




277 






235 






75 






44 




GM Aur 


a 


2xl0~ b 


0.02 


0.3 


2xl0~ b 


0.005 


0.03 


lxl0~ 4 


0.003 


0.02 


4xl0~ 4 


0.002 


0.01 




jcd 


19.2 


20.2 


20.9 


19.2 


20.2 


20.9 


19.2 


19.6 


20.1 


19.2 


19.6 


20.1 




Rc 


4.1 


1000 


4900 


4.1 


915 


4900 


4.1 


45 


233 


4.1 


36 


199 




Ucd 


0.04 


12 


125 


0.04 


11 


125 


0.04 


3.7 


16 


0.04 


3.6 


12 




T c d 


2 


16 


30 


2 


16 


30 


2 


17 


30 


2 


17 


30 




Mcd 


0.55 


1.0 


1.5 


0.55 


1.0 


1.5 


0.55 


0.8 


1.3 


0.55 


0.8 


1.1 




Jcd 


52.3 


53.5 


54.4 


52.3 


53.5 


54.4 


52.3 


52.9 


53.4 


52.3 


52.8 


53.3 




AT 


0.07 


0.6 


4.9 


0.07 


0.6 


4.7 


0.07 


0.4 


3.7 


0.07 


0.5 


3.7 




Age 


1.0 


4.5 


10.0 


1.0 


5.5 


10.0 


1.0 


3.4 


9.8 


1.0 


3.2 


8.6 




AAge 


0.05 


4.8 


9.0 


0.05 


4.3 


8.9 


0.05 


2.4 


8.8 


0.05 


1.5 


5.4 




Mdisfe 


0.03 


0.30 


0.90 


0.03 


0.25 


0.80 


0.03 


0.11 


0.34 


0.03 


0.10 


0.28 




N models 




785 






563 






159 






88 




GM Aur 


/3X10 5 


0.1 


13 


180 


0.1 


4.8 


25 


0.2 


2.2 


10 


0.2 


1.8 


8 




icd 


19.2 


20.2 


20.9 


19.3 


20.2 


20.9 


19.3 


19.8 


20.7 


19.3 


19.8 


20.3 




Rc 


4 


920 


4900 


5 


930 


4900 


5 


105 


1900 


5 


85 


500 




Ucd 


0.02 


10 


120 


0.02 


10 


120 


0.02 


4.5 


30 


0.1 


4.3 


20 




T c d 


2 


15 


30 


2 


16 


30 


2 


15 


30 


2 


15 


30 






0.55 


1.0 


1.5 


0.55 


1.0 


1.5 


0.55 


0.9 


1.3 


0.55 


0.9 


1.2 




Jcd 


52.4 


53.5 


54.4 


52.5 


53.5 


54.4 


52.5 


53.0 


54.1 


52.5 


52.9 


53.5 




AT 


0.08 


0.7 


5 


0.08 


0.6 


5 


0.08 


0.6 


4.8 


0.08 


0.5 


4.3 




Age 


1 


4.4 


10 


1 


5.1 


10 


1 


4.5 


10 


1 


4.4 


10 




AAge 


0.05 


4.9 


8.9 


0.05 


3.8 


8.9 


0.05 


1.6 


6.1 


0.05 


1.0 


3.9 




y^-disk 


0.023 


0.3 


0.87 


0.023 


0.3 


0.80 


0.023 


0.14 


0.80 


0.028 


0.12 


0.35 




N models 




569 






426 






136 






77 





Table 5. Summary of the Monte-Carlo exploration for different sets of constraints. The minimum, maximum and mean value 
of the successful fitting models are given. Symbols: AT Accretion time or molecular cloud collapse time. Units: j cd in cm 2 s _1 , 
i?c in AU, LJcd in units of 10~ 14 s _1 , T c d in K, M c d and Mdisk in solar masses, J c d in g cm 2 s _1 , and AT, Age and AAge in Myrs. 



very large R c also requires a large viscosity, but this time 
so as to yield a large enough accretion rate. As expected, 
no useful constraint on T cc \ can be derived. Constraints 
on the initial rotation of the disk (j cd or w cd ) are weak 



although models with relatively low values of j c d seem to 
be favored. 

The mean derived value for the model parameters to- 
gether with the minimum and maximum allowed values 
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are given on Table 5. The mean value of a inferred from 
this study is 0.02 for DM Tau and 0.002 for GM Aur, both 
close to the standard value of a=0.01 generally found in 
the literature (Hartmann et al, 1998). In the (3 param- 
eterization of turbulence the mean values obtained are 
9 x 10~ 5 and 2 x 10~ 5 for DM Tau and GM Aur respec- 
tively. The (3 value for GM Aur agrees with the prediction 
that j3 w 2 x 10~ 5 from Richard and Zahn (1999). The 
mean value of (3 for DM Tau is a factor of five larger but 
lower values close to the theoretical prediction are allowed. 
These conclusions apply to set {5}, but remain valid for 
sets {2} to 4}. 

A significant part of all models developed gravitational 
instabilities near the end of the disk-formation era and 
mostly for a relatively short time span of 1 — 2 x 10 5 yrs 
in systems with ages less than 3 x 10 5 yrs. These thus 
have a limited impact on the results. We also tested the 
importance of outflows in the DM Tau a case by arbitarily 
imposing an inner boundary at 10 AU. The extra series 
of 500 models show no statistically significant difference 
with the results presented here, except for models with 
very low values of j c( j . 

5.2. Inferred structures of DM Tau and GM Aur 

The selected sets of models provide information not only 
on the regions where we have observational constraints but 
also on the general structure of the disks around DM Tau 
and GM Aur. Figure 12 shows the totality of surface den- 
sity distributions for the four cases studied and different 
sets of constraints. The set of models {3} and {5}, are rep- 
resented as light- and superimposed dark-shaded regions, 
respectively. Fitting the absolute outer radii of the disks 
would provide an additional powerful constraint. This may 
be difficult to obtain since dust observations lack sensitiv- 
ity and CO is affected by photodissociation and conden- 
sation at low temperatures (Aikawa et al, 1999; Dartois 
et al. 2003). Unfortunately, at present it is not clear what 
exactly the outer radii inferred from the millimeter and op- 
tical observations imply: they could either be interpreted 
as an abrupt or as a steady decrease of the density pro- 
file. Without assuming a sharp outer edge for GM Aur, 
models with unreasonable large disk masses become pos- 
sible fits to the observations (set {3}). On the other hand, 
sharp outer edges of circumstellar disks would be in better 
agreement with the (3 prescription of turbulence, as shown 
by fig. 12. 

Temperature is another important quantity whose 
evaluation has important implications for planet forma- 
tion processes. In the a cases it also affects the intensity 
of the viscous turbulence, v, and the disk evolution. The 
envelopes of all possible radial profiles of temperatures 
for the fitted models are shown in Figure 13. These are 
the thermal profiles of our calculated models at the time 
when they verify the selected constraints. As in the pre- 
vious figure, two constraints, {3} and {5}, are shown as 
light- and dark-shaded regions, respectively. At small ra- 



dial distances, the temperature is mostly due to viscous 
dissipation. At large radial distances, the disk becomes op- 
tically thin and the temperatures are essentially controlled 
by the stellar flux. The (3 models arc significatively cooler 
than the a models because they have less mass in the in- 
ner disk, implying lower optical depth and lower central 
temperatures for a given effective (photosphcric) tempera- 
ture. For the same reasons, GM Aur appears to be warmer 
than DM Tau. 

Turbulent viscosity is the key variable to understand 
the disk evolution. Figure 14 shows the envelope of pos- 
sible v profiles as a function of stellar distances for the 
solutions that fit the observations of DM Tau and GM 
Aur (set {5}). At first glance, the solutions for the a and 
[3 models may appear to be relatively similar, a conse- 
quence of the constraint on the density profile provided 
by the observations. However, marked departures between 
the two profiles in the [10— 1000 AU] region are apparent. 
Globally, v a varies as r 3 / 4 , while vp varies as r 1 / 2 , mean- 
ing than the a models are diffused more efficiently in the 
outer disk. This explains why, at large distances from the 
star, the density profile decreases less sharply in the a 
models than in the (3 models. 

5.3. Global evolution of the disk and gravitational 
instabilities 

Figure 15 shows the envelope of possible star and disk 
mass evolution for the models matching conditions {3} 
and {5}, respectively. Up to 15% of our selected models 
(set {5}) were gravitationally unstable, though generally 
only over a short fraction of their evolution, near the end 
of the disk formation era. These correspond to the models 
for which the mass of the disk is close to or exceeds the 
central star mass. The gravitational instability arises at 
30±10 AU and times 1.0 ± 0.6 x 10 5 yr and typically lasts 
for 10 4 - 10 5 yrs for DM Tau and icss than 2 x 10 5 for 
GM Aur due to the larger amount of mass the disk has to 
process. Because we compare observations and models at 
ages larger than 10 6 yrs, gravitational instabilities should 
not affect signficantly our analysis. In sets {4} and {5}, 
models with a > 0.03 were never gravitationally unstable. 
Models with lower values of a had an imposed a — 0.03 
during gravitational instable periods (see section 3.6.3). 

The problem of the evolution of the mass of accretion 
disks as a function of viscosity and initial angular mo- 
mentum J c( j was also studied by Nakamoto and Nakagawa 
(1995). These authors derived a criterion on the viscosity 
parameter a which guarantees that at any given time, the 
disk to star mass ratio never exceeds 0.1: 

log a > -2.0 + 4.51ogJ cd -3.9(logJ cd ) 2 , (29) 

with J c( j in units of 10 53 g cm 2 s" 1 . Their study is relevant 
for systems of ~ 1 M Q collapsing in <~ 6 x 10 5 yrs. They 
considered the parameter range J c d = [6 x 10 52 — 10 53 ] g 
cm 2 s- 1 and a = [10~ 5 - 0.1]. 

For our simulations of DM Tau a we plot on Figure 16, 
the regions of the parameter space that have lower masses 
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Fig. 12. Radial surface density profile of successful models of DM Tau and GM Aur. The dark-shaded regions correspond to 
the models fitting the constraints {5}. The light-shaded regions represent the additional solutions when the constraints are only 
the CO + dust + Accretion rate with large error bar (set {3}). The error bars are the constraints on £ imposed by {5}. 
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Fig. 13. Ensemble of possible mid-plane disk temperatures for two different sets of observational constraints: {3}, light-shaded 
and {5}, dark-shaded. 
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Fig. 14. Viscosities of successful models of DM Tau and GM Aur. Shadowed regions correspond to models calculated with the 
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Fig. 15. Star mass and disk mass as a function of time for models succefully fitting selected observational contraints. The light- 
and dark-shaded regions show the evolution of the disk mass with constraints {3} and {5}, respectively. The dotted and plain 
lines represent the envelopes of possible star masses with the same constraints. The star symbols mark the mean values of the 
ages and disk masses of models fulfilling set {5}. 



than 0.1 M* at the end of the disk formation period (and 
hence, at the moment of maximum disk mass). Our re- 
sults compare well with those of Nakamoto and Nakagawa, 
though their result is strictly valid for systems collaps- 
ing in 0.6 Myr years and ours extend to systems forming 
in 0.05 to 5 Myr. A comparison with fig. 11 shows that 
most solutions for DM Tau and GM Aur are to the left of 
the critical line, i.e. the disk to star mass ratio has been, 
or in some cases is still, larger than 0.1. This seems to 
contradict the fact that most observed disks appear to 



have values of this ratio smaller than 0.1 (e.g. Beckwith 
et al. 1990; Kitamura et al. 2002). However, this contra- 
diction may only be apparent, given the fact that the 
dust absorption coefficient at millimeter wavelengths may 
have been overestimated (see section 2.3), and/or that the 
time spent with a large disk mass is generally relatively 
short. It is also possible that DM Tau and GM Aur are, 
among protoplanetary disks, exceptions, and that most 
other disks have suffered more from either tidal interac- 
tions with nearby stars and/or photoevaporation from ex- 
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M* at the end of the collapse period. The spatial resolution of 
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ternal, relatively massive stars. Alternatively, our assump- 
tions regarding the uncertain opacity factor may be too 
conservative and allow for unrealistically large values of 
the disk mass. This would imply that stronger constraints 
can be derived from the observations presently available. 

5.4. Results for the collapse of magnetized molecular 
clouds 

In order to assess the importance of the collapse model in 
the final results, we also study an a model of DM Tau in 
which the molecular cloud is assumed to collapse under 
the influence of a strong magnetic field. Figure 17 shows 
the final characteristics of models fitting the observational 
constraints when eqs. (9,10) are used. Table 6 shows a 
summary of our results when fitting the observational con- 
straints {2} to {5}. 

Compared with the hydrodynamical collapse mod- 
els, these are characterized by lower values of the ini- 
tial rotational velocity (6 x 10~ 14 s _1 vs. 23 x 10~ 14 
s _1 ) and similar values of the specific angular momentum 
(Logio(jcd) = 20.3 vs. 20.0). This is because differential 
rotation yields a faster rotation of the inner molecular 
cloud core. Furthermore, the final centrifugal radii of the 
disks are smaller than in the non-magnetic case for a given 
rotational velocity. The only possible scenario for the for- 
mation of DM Tau is an inner disk of intermediate size 
(R c ~ 50 AU) that expands in time while it diffuses. This 
is true even if one considers only the weak constraints {2}. 

The distribution of surface densities for the success- 
fully fitting models is shown on fig. 18. Appart from the 
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Fig. 17. Distributions of models fitting the sets of constraints 
{3}, {4} and {5} for DM Tau a in the case of an initially 
magnetized molecular cloud. See fig. 11 for details. 
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Fig. 18. E surface density distributions at different times for 
the successful models in the magnetic cloud case, fulfilling the 
constraints {3} (light-grey) or {5} (dark-grey). See fig. 12 for 
details. 
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Fig. 19. Star and disk mass evolution for the successful models 
in the magnetic cloud case. See fig. 15 for details. 
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System Param. Model results obtained for the different observational constraints 

{2}: CO + Dust {3}: {2} + M {4}: {3} + Outer Limit {5}: {4} + Accurate M 
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Table 6. Summary of Monte-Carlo exploration in the case of an initial collapse of a magnetized molecular cloud core. For each 
set of constraints described in the text {2}, {3} {4} and {5} the minimum, maximum and mean value of the successful fitting 
models are given. Symbols: AT Accretion time or molecular cloud collapse time. Unities: j cd in cm 2 s _1 , R c in AU, co cd in units 
of 10~ 14 s _1 , T c d in K, M c d and Mdisk in solar masses, J c d in g cm 2 s _1 , and, AT, Age and AAge in Myrs. 



fact that high values of the initial angular moment in the 
cloud core are allowed, the results are similar to those ob- 
tained for the non-magnetic case (fig. 13), confirming that 
a strong, sustained viscous evolution of the system is re- 
quired to fit the observations. In particular, the values of 
a derived are almost identical regardless of the collapse 
scenario. 

The global evolution of the star and disk masses is 
shown on fig. 19. Compared to the non-magnetic case 
(fig. 15), the 10 times faster accretion implies that the disks 
grow very rapidly and are likely to become gravitationally 
unstable early on. 

5.5. Turbulence driven by vertical convection? 

One of the advantages of our Monte-Carlo approach is that 
we can test if some of the different theoretical mechanisms 
able to produce turbulence in the disk can explain the ob- 
servations of extended disks. Here we test if the simplest 
characteristics of thermal convection (Lin and Papaloizou, 
1980) can explain the current state of DM Tau. If thermal 
convection is the dominant source of turbulence in cir- 
cumstellar disks, the instability should be essentially lim- 
ited to the inner optically thick disk with a more or less 
sharp transition in the outer optically thin disk (Rudcn 
and Pollack, 1991). 

It is useful to consider some order of magnitude fig- 
ures. In the outer region the opacities are dominated by ice 
grains. For typical temperatures T <~ 10 K, the Rosseland 
opacity (gas+dust) is ~ 0.02 cm 2 g -1 (Pollack et al. 
1986). The disk is unable to sustain a convective verti- 
cal transport of heat at an optical depth t < 1 imply- 
ing a critical surface density T, crit <~ 2/hr ~ 100gcm~ 2 . 
A circumstellar disk of mass 0.1 M Q could then be con- 
vective up to a maximum limit of ~ 50 AU. If we con- 
sider higher temperatures and disk masses (T = 30 K and 



M d = 0.3 M ) this limit extends to 280 AU. Therefore, 
it would be extremely difficult to produce the extended 
structures of disks like DM Tau and GM Aur by the out- 
ward diffusion of initially smaller disks. If turbulent con- 
vection is the dominant source of viscosity, these systems 
should have been formed essentially as they are observed 
now. 

In order to test this possibility we launched a new set of 
calculations for DM Tau following the same range of values 
as in the previous cases and considering an a constant in 
the inner disk but decreasing in the outer part of the disk 
as a function of the optical depth to a given power: 

a = «0 ( — — ) ; T< Tcrit. (30) 

\ T crit/ 

Here ao characterizes the viscosity in the optically thick 
disk, n is an integer number that measures the transition 
of the turbulent active to non active region and T crit — 
1.8 following Rudcn and Pollack (1991). We launched 900 
models with n = 8 (sharp transition) and 500 with n = 1 
(smooth transition) for DM Tau. 

In the first case, n — 8, we obtained 4 models able 
to simultaneously satisfy the CO and dust constrains (set 
{2}) with centrifugal radius around 1300 AU. These mod- 
els were typically far too massive and could not satisfy 
simultaneously the observed accretion rate. 

In the second case, we chose n = 1. This situation 
corresponds to a slow decrease of turbulence proportional 
to the optical depth. It is probably unrealistic but could 
possibly model a hybrid case in which convection is active 
at r > 1 and another undefined, weaker, turbulent source 
sets in at lower optical depths. In this case, material can be 
transported from the optically thin to the optically thick 
regions. 

We found that 19 % of all models are able to fit the 
CO + dust requirements (set {2}) and that 13 % also 
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Fig. 20. Global evolution of the star and disk mass. This was 
the only case found where turbulence of convective origin was 
found to satisfy the most retrictive set of constraints {5}. See 
fig. 5 for details. 
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Fig. 21. £ surface density distributions at different times for 
the same model as the previous figure. See fig. 6 for details. 



fit the accretion rate constraint (set {3}). This is simi- 
lar to the results obtained previously with a uniform a 
value. However, the models fitting the constraints {3} 
are now characterized by larger values of the centrifu- 
gal radius (R c = [45,4900] AU), specific angular mo- 
mentum (Logio(jcd) = [19.7 — 20.8]) and hence viscosity 
(a = [0.002,0.6]) than those listed in Table 5. These disks 
are generally massive, they develop gravitational instabil- 
ities and are formed basically as they are observed. We 
found that only a tiny fraction ~ 2% of all models also fit 
constraints {4}. These are characterized by: R c — [45, 300] 
AU; a = [0.003,0.07]; Log w (j c d) = [19.7 - 20.1]. Finally, 
only one single model was able to also fit simultaneously 
the strict accretion rate we have considered on set {5}. 

The global evolution of this model together with the 
radial distributions of its surface density and the disk ef- 
fective viscosity are shown in figs. 20, 21 and 22. The 
model is characterized by a = 0.07, j c d = 19.73, R c = 
62 AU and total accretion time 0.22 Myr. The molecu- 
lar cloud parameters are also found to agree reasonably 
well with those inferred for the Taurus formation region 
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Fig. 22. Radial distribution of turbulent viscosity in the disk 
as a function of time. As the disk spreads it is able to develop 
viscosity even in regions far beyond the final centrifugal radius 
R c =62 AU. This is because this case corresponds to n = 1 (see. 
eq. (30)), implying that the transition between the convective 
to non-convective regime is extremely smooth (and probably 
unrealistically so). 

(uj cd = 4.38 x 10~ 14 s" 1 , T cd = 12 K, M = 0.05 M and 
M c d = 0.527 M Q ). It must be stressed however that this 
chance result requires initial conditions that are very close 
to those listed here. 

We conclude that it is extremely unlikely that convec- 
tive instabilities could have been the dominant source of 
angular momentum transport in DM Tau and GM Aur. 

6. Conclusions 

We presented a (relatively) simple model of evolution of 
circumstellar disks and applied this model systematically 
to two well-characterized objects: DM Tau and GM Aur. 
We chose to survey extensively a rather large parameter 
space, using various observational constraints from CO 
observations of these disks, dust emission properties, de- 
rived accretion rates and age and masses of the central 
stars. Two viscosity parameterizations, a and (3, and two 
cloud collapse models were used. 

We showed that the current state of DM Tau and GM 
Aur can be understood as resulting from the collapse of a 
molecular cloud core, the formation of an accretion disk 
and its spreading due to angular momentum conservation 
and turbulent diffusive transport. This scenario is con- 
sistent with a specific angular momentum initially in the 
cloud core (j c d ~ 10 20 — 10 21 cm 2 s _1 ) similar to measure- 
ments in class I protostellar cores (Ohashi et al. 1997), and 
values of the viscosity that are within an order of magni- 
tude of the standard a ~ 0.01 derived from measured ages 
and accretion rates (e.g. Hartmann et al. 1998). The signif- 
icant outward diffusion implies that most of the material 
has been reprocessed and lies far from the position where 
it had originally fallen onto the disk. 

More specifically, the values of a that can be inferred 
from this study range within 0.001 — 0.1 for DM Tau and 
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4 x 1(T 4 - 0.01 for GM Aur. In the case of the (3 models, 
the likely values of this parameter are 2 x 10~ 5 — 5 x 10 -4 
for DM Tau and 2 x 10~ 6 - 8 x 10~ 5 for GM Aur. 

Unfortunately, the still relatively large error bars on 
observationally determined quantities (in particular the 
age of the system and the star accretion rate) prevent 
inferring tight constraints on parameters characterizing 
the systems. We cannot argue in favor of either the a 
or /? parameterizations of turbulent viscosity. Similarly, 
our results are consistent with both a relatively slow (hy- 
drodynamic) or a fast (magnetic) collapse of the primor- 
dial molecular cloud core. We can however rule out ther- 
mal convection as the main source of turbulence in the 
disk: The observed disks of DM Tau and GM Aur are 
too large and extend too far in the optically thin regime 
in which convection would be suppressed. On the other 
hand, they cannot be formed directly from the molecu- 
lar cloud. Convection could still play a role in the inner 
disk, but another source of turbulence is required in the 
optically thin regime. 

The current state of both systems could be explained 
by a wide variety of evolutionary histories. In most cases, 
the disks were always gravitationally stable. However, for 
low- and high-values of the initial angular momentum j c( j 
(corresponding to R c < 10 AU or R c > 2000 AU), disks 
can become gravitationally unstable for <~ 10 5 yrs or less. 
In the case of DM Tau and GM Aur, gravitational insta- 
bilities thus appear to have had a limited role. 

The relatively large values of the viscosity that we de- 
rive both for DM Tau and GM Aur imply that an efficient 
turbulent diffusion mechanism is present throughout these 
disks. As a consequence, any giant planet forming in these 
disks will be in risk of migrating rapidly into its star, ex- 
cept if (i) it is protected by an inner dead zone character- 
ized by a low viscosity or (ii) it forms late, i.e. when the 
local disk mass is comparable to the mass of the planet 
itself. We regard the second possibility as more likely, but 
this requires more specific studies. In any case, the pres- 
ence of extended viscous disks around DM Tau and GM 
Aur has consequences on how we should apprehend planet 
formation and should be included in future studies of this 
problem. 

Several improvements to this work are possible. 
Obviously, our sets of constraints {1} to {5} should even- 
tually be replaced by a global assessment of the compat- 
ibility of the derived density profiles with the different 
observations (SED, dust cmissivity, CO observations, re- 
flected light observations). This is not trivial because of 
possible variations of the dust to gas ratio (dust migration 
due to gas drag, condensation at different temperatures), 
of radial and vertical variations of the size distribution 
of solid particles, and of variations in the chemical com- 
position (in particular due to CO condensation and its 
photodissociation). Another improvement would be to in- 
clude the thermal evaporation of the disks at large orbital 
distances (hundreds of AU). 

In any case, we believe that meaningful constraints on 
physical processes leading to the formation of extended 



disks can now be derived. These critically depend on the 
inferred surface density profiles at 100 AU and beyond. 
High sensitivity observations able to detect small amounts 
of gas or dust in the outer regions of the disks such as 
forseen for ALMA are of course highly desirable. A better 
understanding of the inner disk regions and constraints on 
the rate of accretion onto the star would also be very valu- 
able. Ongoing programs (in particular observations with 
the Spitzer telescope, but also from ground-based facili- 
ties) should ensure a rapid progress towards understand- 
ing the formation and evolution of protoplanetary disks. 
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